/*************************************************************************
Purpose:  This program loads the input file after it has been
          processed by stack_read_initial.py, then loads the orbit files
          from Yuan.  Next, the orbit files are transformed into 
          geocentric (or geodetic?) coordinates, and the datapoints are
          synced with the input file.  A new input file is produced which
          contains both sets of data.
**************************************************************************
Input: The input file and the orbit files are read from disk.
**************************************************************************
Output: A new input file with orbit data included is output.
**************************************************************************
Copyright 2010 Bryan Killett

This file is part of GAIA.

GAIA is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

GAIA is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with GAIA. If not, see <http://www.gnu.org/licenses/>.
*************************************************************************/
//Header file "definitions.hpp" declares many functions, constants, datatypes.
//In order to use the functions declared in definitions.hpp, must link
//with object code that includes functions.cpp.
#include "definitions.hpp"

int main(){

  input_s input;
  long long i,j;
  char s[max_length];
  string basefolder(gaiafolder());
  string file;
  FILE *out_fp;
  long long time_loop,time_total;//Used for estimated time left.
  time_t start_time,current_time=-1,previous_time=-1;//Used to size notif.   
  long long orbit_points;//Holds the number of orbit points in the current file.
  long long current_orbit_counter;//Runs from 0 to orbit_points-1.

  //Since this program is run without every loading the entire input file
  //(which is intentional!) I need to set a variable containing the estimated
  //number of input points in order to make the estimate of time remaining
  //correct.
  //long long input_numpoints_estimated = 6500000;
  long long input_numpoints_estimated = 9000000;
  double gap=15;//Minimum time b/w orbit pts to qualify as a data gap.
  input.options.segmentsize = 100000;//Load "segmentsize" input pts at a time.
  
  //Used for converting timespans in seconds to human readable units.
  timespan_conversion_s timespan_conversion;
  
  start_time = time(NULL);//Record the time at the start of the program.

  ////////////////////////////////////////////////////////////////////////
  //Load first segment of input data from disk.
  ////////////////////////////////////////////////////////////////////////  
  input.options.basefolder= "data/stack_output/deg150/";
  //Input accelerations, midpoint lats/lons.
  input.options.filename = "stack_deg150_2008_20s_smoothing_50N.txt";
  //Output filename of the final input and orbit combined file.
  //input.options.output_filename = "stack_and_orbit_7yrs_new.txt";
  input.options.output_filename="stack_and_orbit_deg150_2008_20s_smoothing_50N.txt";
  
  long long mainindex = 0;//We're starting to read the main input at index "0".
  //Need an initialize parameter.
  //This parameter can also double as a "main" input index - doesn't get
  //reset when the next chunk of input data is read again.  It marks the
  //index in the MAIN input file of the next input point to be read
  input = load_another_input_segment(input.options,mainindex);
  long long segmentcounter = 1;//Already loaded 1 input segment.
  //Save file marking progress (useful when run remotely).
  file = gaiafolder();
  file.append(input.options.basefolder);
  sprintf(s,"combined_%06lld_segments",segmentcounter-1);
  file.append(s);
  out_fp = fopen(file.c_str(),"w");
  fclose(out_fp);

  ////////////////////////////////////////////////////////////////////////
  //Load orbit files one by one, while syncing with the input.
  //Only transform coordinates when the orbit time syncs with input time.
  ////////////////////////////////////////////////////////////////////////  
  
  //These variables are used to search the orbit files folder.
  DIR *orbitDIR;
  struct dirent *orbitDirEnt;  
  string orbitfolder(gaiafolder());
  DIR *orbitsubDIR;
  struct dirent *orbitsubDirEnt;  
  string orbitsubfolder;
  char month[3];
  int numfiles,day,temp;
  //This vector is used to sort the subfolders within the orbit folder.
  vector<string>subfolder_names;
  int subfolder_index;//Used to step through sorted subfolder vector.
    
  //These variables are used to open orbit files.
  vector<int> date,date_original;
  vector<string> filename,filename_original;
  
  //These variables are used to sync the orbit files with the input.
  double st;//Current input time.
  double st_minus_ot;//Current input time minus current orbit time.  
  double dt_current;//Absolute value of st_mins_ot - delta time (current).
  long long si;//Current input index.
  int dummy = -1,found;
  string file2;//Holds the filename of a grace2 orbit file.
  string temp_string, temp_string1, temp_string2;
  double t1,x1,y1,z1,t2,x2,y2,z2;
  vector<double> file_t1, file_x1, file_y1, file_z1, file_x2, file_y2, file_z2;
  double vx,vy,vz;
  double xp1,yp1,zp1,xp2,yp2,zp2;//interpolated values for x,y,z for grace1,2.
  double xm,ym,zm;//Midpoint between the satellites, based on interp. values  
  double latm,lonm,altm;//Midpoint in geodetic (centric) coords.
  double dt[3],orbit_times[3],x1_array[3],y1_array[3],z1_array[3],x2_array[3],y2_array[3],z2_array[3];
  long long orbit_file_counter;//Counts the number of orbit files that have been opened.
  long long closest_points_array_filled;//Set to 1 if array is filled, 0 to start.
  int print_debug = 0;//Controls various debug messages.  Set to 0 to silence.
  int skip_first_file=1;//Used to skip 1st orbit file b/c it's just pre-grace.
  //Initialize variables.
  si = 0;
  st = (double)input.time.at(si);
  if(print_debug>0) cout <<"Orbit times closest to input.time["<<si<<"]= "<<input.time.at(si)<<" are: ";
  dt[0] = 100000000;
  dt[1] = 200000000;
  dt[2] = 300000000;
  orbit_times[0] = -1;
  orbit_times[1] = -2;
  orbit_times[2] = -3;
  orbit_file_counter = 0;    
  closest_points_array_filled = 0;
    
  //Open the orbit folder.
  orbitfolder.append("data/orbit_files");
  orbitDIR = opendir(orbitfolder.c_str());

  //Examine the main orbit folder for subfolders such as 2002,2003, etc.
  //Read all of them at once, then sort them first to insure that they're
  //read in the order 2002,2003,etc.
  orbitDirEnt = readdir(orbitDIR);
  while (orbitDirEnt != NULL){
    //Open the subfolder inside the main orbit folder.
    orbitsubfolder=orbitfolder;
    orbitsubfolder.append("/");
    strcpy(s,orbitDirEnt->d_name);
    orbitsubfolder.append(s);
    //If this entry isn't irrelevant, add it to the subfolder_names vector.
    if(strcmp(s,".") and strcmp(s,"..") and strcmp(s,"orbit_notes.txt") and strcmp(s,"orbit_notes.txt~") and strcmp(s,".gz") and strcmp(s,"temp") and strcmp(s,"old")){ 
      subfolder_names.push_back(orbitsubfolder);
    }
    orbitDirEnt = readdir(orbitDIR);    
  }
  
  //Close the main orbit folder.
  closedir(orbitDIR);  
  
  //Now that all the orbit subfolder names have been read, sort them.
  sort(subfolder_names.begin(),subfolder_names.end());

  //Debug.
  for(i=0;i<(long long)subfolder_names.size();i++){
    cout<<"After sorting, subfolder #"<<i+1<<" of "<<subfolder_names.size()<<": "<<subfolder_names.at(i)<<endl;
  }    
      
  //Once subfolders have been sorted, open each one and read the orbit files
  //contained within.
  for(subfolder_index=0;subfolder_index<(int)subfolder_names.size();subfolder_index++){
    //Open the subfolder inside the main orbit folder.
    orbitsubfolder = subfolder_names.at(subfolder_index);
    cout<<"Working in folder "<<orbitsubfolder<<endl;    
    //Set the number of files in this subfolder to 0 and clear vectors.
    numfiles = 0;
    date.clear();
    filename.clear();
    orbitsubDIR = opendir(orbitsubfolder.c_str());
    //Examine the orbit subfolder and record the filenames, along
    //with creating a vector that orders the files in chronological order.
    orbitsubDirEnt = readdir(orbitsubDIR);
    while (orbitsubDirEnt != NULL){
      strcpy(s,orbitsubDirEnt->d_name);       
      //Only open grace1 files- no need to open grace2 files yet.
      if(strcmp(s,".") and strcmp(s,"..") and strstr(s,"grace2") == NULL){         
        //Extract the date from the filename, convert it to a 
        //number which (sorta) represents the date, but always sorts right.
        sscanf(s,"grace1_p_%02d-%3s", &day,&month);
        if(!strcmp(month,"JAN")){
          temp = 100;
        }
        if(!strcmp(month,"FEB")){
          temp = 200;
        }
        if(!strcmp(month,"MAR")){
          temp = 300;
        }
        if(!strcmp(month,"APR")){
          temp = 400;
        }
        if(!strcmp(month,"MAY")){
          temp = 500;
        }
        if(!strcmp(month,"JUN")){
          temp = 600;
        }
        if(!strcmp(month,"JUL")){
          temp = 700;
        }
        if(!strcmp(month,"AUG")){
          temp = 800;
        }
        if(!strcmp(month,"SEP")){
          temp = 900;
        }
        if(!strcmp(month,"OCT")){
          temp = 1000;
        }
        if(!strcmp(month,"NOV")){
          temp = 1100;
        }
        if(!strcmp(month,"DEC")){
          temp = 1200;
        }
        temp += day;
        date.push_back(temp);
        filename.push_back(s);
        numfiles += 1;
      }
      //Read the next filename for the next loop.
      orbitsubDirEnt = readdir(orbitsubDIR);      
    }
    date_original = date;
    filename_original = filename;

    //Now, sort the dates of files in this folder.
    sort(date.begin(),date.end());
    
    //Use the sorted dates to sort the filenames, caveman style!
    for(i=0;i<numfiles;i++){
       j=0;
       found=0;
       while(found==0 and j<numfiles){
         if(date_original.at(j) == date.at(i)) found = 1;
         else j++;
       }
       filename.at(i) = filename_original.at(j);
    }
      
    //Read data from each file in chronological order.
    if(skip_first_file > 0) cout <<"Skipping "<<skip_first_file<<" file(s) at beginning of first year only."<<endl;
    for(i=skip_first_file;i<numfiles;i++){//Set i=1 in 1st folder to skip 1ST orbit file.
      printf("Loading data from file %lld of %d\n",i+1,numfiles);
      orbit_file_counter++;
      //Load grace 1 and grace2 files for the same time period.
      temp_string = orbitsubfolder + "/"+filename.at(i);
      ifstream grace1_file (temp_string.c_str());
      file2 = temp_string;
      file2.replace(file2.find("grace1"),6,"grace2");
      ifstream grace2_file (file2.c_str());
      //Before loading file, clear vectors.
      file_t1.clear();
      file_x1.clear();
      file_y1.clear();
      file_z1.clear();
      file_x2.clear();
      file_y2.clear();
      file_z2.clear();
      while(getline(grace1_file,temp_string1) and getline(grace2_file,temp_string2)){
        //Examine each line of each file and extract time and positions.
        //First, change scientific notation character from "D" to "e".
        for(j=0;j<7;j++){
          temp_string1.replace(temp_string1.find("D"),1,"e");
          temp_string2.replace(temp_string2.find("D"),1,"e");
        }
        sscanf(temp_string1.c_str(),"%le %le %le %le %le %le %le",&t1,&x1,&y1,&z1,&vx,&vy,&vz);
        sscanf(temp_string2.c_str(),"%le %le %le %le",&t2,&x2,&y2,&z2);
        //Transform time to seconds past J2000.
        t1 = (t1 - 2451545.0)*86400.0;//Subtract JD of J2000, then from days to seconds.
        //Save to vectors.
        file_t1.push_back(t1);
        file_x1.push_back(x1);
        file_y1.push_back(y1);
        file_z1.push_back(z1);
        file_x2.push_back(x2);
        file_y2.push_back(y2);
        file_z2.push_back(z2);
      }
      orbit_points = file_t1.size();
      grace1_file.close();
      grace2_file.close();
      cout<<"Finished loading "<<orbit_points<<" orbit points from file "<<filename.at(i)<< endl;

      if(dummy == -1) dummy=0;//Holds # o orbit times since last gap.
      for(current_orbit_counter=0; current_orbit_counter<orbit_points; current_orbit_counter++){
        if(dummy == 3) {
          closest_points_array_filled = 1;
        }
        dummy++;//Need to fill closest points array first!  Used after gaps.

        //Load orbit data from vectors for this loop.
        t1 = file_t1.at(current_orbit_counter);
        x1 = file_x1.at(current_orbit_counter);
        y1 = file_y1.at(current_orbit_counter);
        z1 = file_z1.at(current_orbit_counter);
        x2 = file_x2.at(current_orbit_counter);
        y2 = file_y2.at(current_orbit_counter);
        z2 = file_z2.at(current_orbit_counter);
        
        //Compare this time to the time in input index si.
        start_of_comparison://Label allowing goto statement to return here.
        
        st_minus_ot = st - t1;
        //Calculate absolute value, store in dt_current.
        dt_current = fabs(st_minus_ot);
        //if(1==2){
        //  printf("input.time[%08lld]: %13.5f\n",si,st);
        //  printf("orbit time          : %13.5f. Also, st - ot:%13.5f\n",t1,st-t1);
        //  printf("dt[0]: %13.5f\n",dt[0]);
        //  printf("dt[1]: %13.5f\n",dt[1]);
        //  printf("dt[2]: %13.5f\n",dt[2]);
        //  printf("orbit_times[0]: %13.5f\n",orbit_times[0]);
        //  printf("orbit_times[1]: %13.5f\n",orbit_times[1]);
        //  printf("orbit_times[2]: %13.5f\n",orbit_times[2]);
        //  cout<<endl;
        //}
        //Is this the closest orbit point (relative to "st") ever?
        if(dt_current < dt[0]){
          //cout<<"This is the closest point ever!"<<endl;
          //If so, swap 2nd and 3rd closest values, 
          //then swap 1st and 2nd values, then store current in 1st.
          for(j=2;j>0;j--){
            swap_double(&dt[j],&dt[j-1]);
            swap_double(&orbit_times[j],&orbit_times[j-1]);
            swap_double(&x1_array[j],&x1_array[j-1]);
            swap_double(&y1_array[j],&y1_array[j-1]);
            swap_double(&z1_array[j],&z1_array[j-1]);
            swap_double(&x2_array[j],&x2_array[j-1]);
            swap_double(&y2_array[j],&y2_array[j-1]);
            swap_double(&z2_array[j],&z2_array[j-1]);
          }
          dt[0] = dt_current;
          orbit_times[0] = t1;
          x1_array[0] = x1;
          y1_array[0] = y1;
          z1_array[0] = z1;
          x2_array[0] = x2;
          y2_array[0] = y2;
          z2_array[0] = z2;
          //printf("input.time[%08lld]: %13.5f\n",si,st);
          //printf("orbit time          : %13.5f. Also, st - ot:%13.5f\n",t1,st-t1);
          //printf("dt[0]: %13.5f\n",dt[0]);
          //printf("dt[1]: %13.5f\n",dt[1]);
          //printf("dt[2]: %13.5f\n",dt[2]);
          //printf("orbit_times[0]: %13.5f\n",orbit_times[0]);
          //printf("orbit_times[1]: %13.5f\n",orbit_times[1]);
          //printf("orbit_times[2]: %13.5f\n",orbit_times[2]);
          //cout<<endl;
          
        }//End of "1st closest orbit point ever" if statement.
        //If not, is this the 2nd closest point?
        else if(dt_current < dt[1]){
          //cout<<"This is the second closest point ever!"<<endl;
          //If so, swap 2nd and 3rd closest values, put current
          //values in 2nd slot.  
          for(j=2;j>1;j--){
            swap_double(&dt[j],&dt[j-1]);
            swap_double(&orbit_times[j],&orbit_times[j-1]);
            swap_double(&x1_array[j],&x1_array[j-1]);
            swap_double(&y1_array[j],&y1_array[j-1]);
            swap_double(&z1_array[j],&z1_array[j-1]);
            swap_double(&x2_array[j],&x2_array[j-1]);
            swap_double(&y2_array[j],&y2_array[j-1]);
            swap_double(&z2_array[j],&z2_array[j-1]);              
          }  
          dt[1] = dt_current;
          orbit_times[1] = t1;
          x1_array[1] = x1;
          y1_array[1] = y1;
          z1_array[1] = z1;
          x2_array[1] = x2;
          y2_array[1] = y2;
          z2_array[1] = z2;
          //printf("input.time[%08lld]: %13.5f\n",si,st);
          //printf("orbit time          : %13.5f. Also, st - ot:%13.5f\n",t1,st-t1);
          //printf("dt[0]: %13.5f\n",dt[0]);
          //printf("dt[1]: %13.5f\n",dt[1]);
          //printf("dt[2]: %13.5f\n",dt[2]);
          //printf("orbit_times[0]: %13.5f\n",orbit_times[0]);
          //printf("orbit_times[1]: %13.5f\n",orbit_times[1]);
          //printf("orbit_times[2]: %13.5f\n",orbit_times[2]);
          //cout<<endl;
          
        }//End of "2nd closest orbit point" if statement.
        //If the current orbit time is greater than the input time,
        //and the closest points array is filled (i.e. at least 3
        //loops have been performed, and the current orbit time is NOT
        //the 1st or 2nd closest, that means we've already got the closest
        //2 times stored in the array.  Interpolate, push_back() coords and
        //then advance to the next input time.  Do this only if the first
        //two delta times are less than "gap".
        else if(st_minus_ot < 0 and closest_points_array_filled == 1 and dt_current > dt[0] and dt_current > dt[1] and dt[0] < gap and dt[1] < 2*gap){

          //cout<<"These orbit times are close enough to be interpolated!"<<endl;          
          //printf("input.time[%08lld]: %13.5f\n",si,st);
          //printf("orbit time          : %13.5f. Also, st - ot:%13.5f\n",t1,st-t1);
          //printf("dt[0]: %13.5f\n",dt[0]);
          //printf("dt[1]: %13.5f\n",dt[1]);
          //printf("dt[2]: %13.5f\n",dt[2]);
          //printf("orbit_times[0]: %13.5f\n",orbit_times[0]);
          //printf("orbit_times[1]: %13.5f\n",orbit_times[1]);
          //printf("orbit_times[2]: %13.5f\n",orbit_times[2]);
          //cout<<endl;
        
          if(print_debug>0) cout <<orbit_times[0]<<" and "<<orbit_times[1]<<endl;
          //Whenever a new point isn't closer than the previously closest point, the program is about to overshoot the closest two points.
          //Interpolate based on these times and distances.
          double time_factor = (st-orbit_times[1])/(orbit_times[1]-orbit_times[0]);
          xp1 = (x1_array[1]-x1_array[0])*time_factor + x1_array[1];
          xp2 = (x2_array[1]-x2_array[0])*time_factor + x2_array[1];
          yp1 = (y1_array[1]-y1_array[0])*time_factor + y1_array[1];
          yp2 = (y2_array[1]-y2_array[0])*time_factor + y2_array[1];
          zp1 = (z1_array[1]-z1_array[0])*time_factor + z1_array[1];
          zp2 = (z2_array[1]-z2_array[0])*time_factor + z2_array[1];                                                            
          xm = (xp1+xp2)/2; ym = (yp1+yp2)/2; zm = (zp1+zp2)/2;            
          
          //Transform midpoint to lat,lon,alt to compare to input.
          xyz2llh(xm,ym,zm,latm,lonm,altm);
          double km_error = 110*180/pi*acos(cos((90-input.lat.at(si))*pi/180)*cos((90-latm)*pi/180) + sin((90-input.lat.at(si))*pi/180)*sin((90-latm)*pi/180)*cos((input.lon.at(si)-lonm)*pi/180));
          if(km_error > 0.35){
            cout<<"!!!HUGE Midpoint error (km): "<<km_error<<" (new: "<<spherical_distance(input.lat.at(si),input.lon.at(si),latm,lonm)*planet_avg_radius<<") at time "<<input.time.at(si)<<", lat "<<input.lat.at(si)<<", and lon "<<input.lon.at(si)<<endl;
            //Mark this input time for deletion.
            input.bad_indices.push_back(si);
          }

          //Save interpolated positions as values in the input structure.
          input.xA.push_back(xp1);
          input.xB.push_back(xp2);
          input.yA.push_back(yp1);
          input.yB.push_back(yp2);
          input.zA.push_back(zp1);
          input.zB.push_back(zp2);
          
          //Now go to the next input time.
          si++;
          
          //If si indicates that the end of input.time has been reached,
          //either load more data from the input file or, if it's empty,
          //jump out of the loop by writing a (gasp!) safe goto function
          //that first closes the open files.
          if(si == (long long)input.time.size()){
            if(mainindex>0){//i.e. if more input points remain to be read.
            
              //Calculate time since start of program and print
              //an estimate of the time remaining.
              if(current_time != -1) previous_time = current_time;
              else previous_time = start_time;
              current_time = time(NULL);
              time_loop = current_time - previous_time;
              time_total = current_time - start_time;
              timespan_conversion.seconds=time_total;
              printf("%5.2f%% finished if st.np==%7lld.",(double)(segmentcounter*input.options.segmentsize)*100.0/(double)input_numpoints_estimated,input_numpoints_estimated);
              cout<<endl;

              if(print_debug > 0) cout<<"Mainindex: "<<mainindex<<endl;
              append_input_and_orbit_data(input,mainindex);
              input = load_another_input_segment(input.options,mainindex);
              si = 0;
              //Save file marking progress (useful when run remotely).
              file = gaiafolder();
              file.append(input.options.basefolder);
              sprintf(s,"combined_%06lld_segments",segmentcounter);
              file.append(s);
              out_fp = fopen(file.c_str(),"w");
              fclose(out_fp);
              segmentcounter++;
            }
            else{  
              append_input_and_orbit_data(input,mainindex);
              grace1_file.close();
              grace2_file.close();
              cout<<"End of input file reached.\n"<<endl;
              goto end_of_orbit_files;
            }
          }//End of "if end of input structure has been reached" statement.
          st = (double)input.time.at(si);

          //Recalculate dt's to the three times stored in time[0-2] based on new input time.
          dt[0] = fabs(st - orbit_times[0]);
          dt[1] = fabs(st - orbit_times[1]);
          dt[2] = fabs(st - orbit_times[2]);
          
          //Reorder "closest" arrays based on new dts.
          for(j=2;j>0;j--){
            if(dt[j] < dt[j-1]){
              swap_double(&dt[j],&dt[j-1]);
              swap_double(&orbit_times[j],&orbit_times[j-1]);
              swap_double(&x1_array[j],&x1_array[j-1]);
              swap_double(&y1_array[j],&y1_array[j-1]);
              swap_double(&z1_array[j],&z1_array[j-1]);
              swap_double(&x2_array[j],&x2_array[j-1]);
              swap_double(&y2_array[j],&y2_array[j-1]);
              swap_double(&z2_array[j],&z2_array[j-1]);
            }
          }
          //Jump back to previous point in loop right after loading the
          //next orbit time to recalculate everything for current input
          //time as well.
          goto start_of_comparison;
          
        }//End of "time to interpolate and save" if statement.
        //If none of these cases are true, something strange is
        //happening, like a gap in the orbit files.
        else{
          if(dt_current < dt[2]){
            //If so, put current value in 3rd slot.  
            //cout<<"Current point is third closest."<<endl;
            //printf("input.time[%08lld]: %13.5f\n",si,st);
            //printf("dt[0]: %13.5f\n",dt[0]);
            //printf("dt[1]: %13.5f\n",dt[1]);
            //printf("dt[2]: %13.5f\n",dt[2]);
            //printf("orbit_times[0]: %13.5f\n",orbit_times[0]);
            //printf("orbit_times[1]: %13.5f\n",orbit_times[1]);
            //printf("orbit_times[2]: %13.5f",orbit_times[2]);
            //cout<<endl;

            //cout<<"orbit_file_counter: "<<orbit_file_counter<<endl;
            //cout<<"input.time index: "<<si<<endl;
            dt[2] = dt_current;
            orbit_times[2] = t1;
            x1_array[2] = x1;
            y1_array[2] = y1;
            z1_array[2] = z1;
            x2_array[2] = x2;
            y2_array[2] = y2;
            z2_array[2] = z2; 
          }
          
          //Is this input point in an orbit data gap?
          //Check to see if the delta times are uniformly larger than the
          //"gap".  If so, increase the input time and ERASE the input
          //time that doesn't have orbit data!
          if (dt[0] > gap and dt[1] > gap and dt[2] > gap){
            //cout<<"WARNING: input time is in an orbit data gap!"<<endl;
            //cout<<"orbit_file_counter: "<<orbit_file_counter<<endl;
            //printf("soon to be deleted input.time[%08lld]: %13.5f\n",si,st);
            //printf("dt[0]: %13.5f\n",dt[0]);
            //printf("dt[1]: %13.5f\n",dt[1]);
            //printf("dt[2]: %13.5f\n",dt[2]);
            //printf("orbit_times[0]: %13.5f\n",orbit_times[0]);
            //printf("orbit_times[1]: %13.5f\n",orbit_times[1]);
            //printf("orbit_times[2]: %13.5f\n",orbit_times[2]);
            //cout<<"Deleting input elements with index "<<si<<"."<<endl;
            cout<<"Deleting input pt with time "<<input.time.at(si)<<", lat "<<input.lat.at(si)<<", and lon "<<input.lon.at(si)<<endl;//Subtract
            input.time.erase(input.time.begin()+si);
            input.lat.erase(input.lat.begin()+si);
            input.lon.erase(input.lon.begin()+si);
            input.accel.erase(input.accel.begin()+si);
            input.ascending.erase(input.ascending.begin()+si);
            
            //Do NOT increment si- now it points to the next input time
            //because the si'th elements just got erased.
            //If si indicates that the end of input.time has been reached,
            //either load more data from the input file or, if it's empty,
            //jump out of the loop by writing a (gasp!) safe goto function
            //that first closes the open files.
            if(si == (long long)input.time.size()){
              if(mainindex >0){
              
                //Calculate time since start of program and print
                //an estimate of the time remaining.
                if(current_time != -1) previous_time = current_time;
                else previous_time = start_time;
                current_time = time(NULL);
                time_loop = current_time - previous_time;
                time_total = current_time - start_time;
                timespan_conversion.seconds=time_total;
                printf("%5.2f%% finished if st.np==%7lld.  Est. time left:",(double)(segmentcounter*input.options.segmentsize)*100.0/(double)input_numpoints_estimated,input_numpoints_estimated);
                cout<<sec2human(timespan_conversion)<<endl;
              
                if(print_debug > 0) cout<<"Mainindex: "<<mainindex<<endl;
                append_input_and_orbit_data(input,mainindex);
                input = load_another_input_segment(input.options,mainindex);
                si = 0;
                //Save file marking progress (useful when run remotely).
                file = gaiafolder();
                file.append(input.options.basefolder);
                sprintf(s,"combined_%06lld_segments",segmentcounter);
                file.append(s);
                out_fp = fopen(file.c_str(),"w");
                fclose(out_fp);  
                segmentcounter++;
              }
              else{  
                append_input_and_orbit_data(input,mainindex);
                grace1_file.close();
                grace2_file.close();
                cout<<"End of input file reached.\n"<<endl;
                goto end_of_orbit_files;
              }
            }//End of "if end of input struct has been reached" if statement
            st = (double)input.time.at(si);
            //printf("new input.time[%08lld]: %13.5f",si,st);
            cout<<endl;              
            
            //Recalculate dt's to the three orbit times stored in orbit_times[0-2] based on new input time.
            dt[0] = fabs(st - orbit_times[0]);
            dt[1] = fabs(st - orbit_times[1]);
            dt[2] = fabs(st - orbit_times[2]);
            
            //Reorder "closest" arrays based on new dts.
            for(j=2;j>0;j--){
              if(dt[j] < dt[j-1]){
                swap_double(&dt[j],&dt[j-1]);
                swap_double(&orbit_times[j],&orbit_times[j-1]);
                swap_double(&x1_array[j],&x1_array[j-1]);
                swap_double(&y1_array[j],&y1_array[j-1]);
                swap_double(&z1_array[j],&z1_array[j-1]);
                swap_double(&x2_array[j],&x2_array[j-1]);
                swap_double(&y2_array[j],&y2_array[j-1]);
                swap_double(&z2_array[j],&z2_array[j-1]);
              }
            }

            //Jump back to previous point in loop right after loading the
            //next orbit time to recalculate everything for current input
            //time as well.
            goto start_of_comparison;
          }//End of "is this point in a data gap" if statement.
          
        }//End of "if orbit pts aren't 1st,2nd or close enough for interp.
      }//End of loop through data from the current file.
    }//End of loop through files in this subfolder.
    skip_first_file = 0;//Don't skip first file other than in 2002 folder.
    //Close the orbit subfolder.
    closedir(orbitsubDIR);
    //Clear vectors used in this folder.
    date_original.clear();
    date.clear();
    filename.clear();
  }//End of loop through subfolders.
  
  //When the end of the input file is reached, the program goes here:
  end_of_orbit_files:
  
  //Have the number of input points changed?  Hopefully not, but check!
  if(1==2){//Now that input loads in segments, this is meaningless.
    if((long long)input.time.size() != (long long)input.xA.size()){
      cout<<"numpoints != xA.size()!"<<endl;
      cout<<"input.time.size() = "<<input.time.size()<<" and input.xA.size() = "<<input.xA.size()<<endl;
    }  
    if((long long)input.time.size() != (long long)input.time.size()){
      cout<<"numpoints != time.size()!"<<endl;
      cout<<"input.time.size() = "<<input.time.size()<<" and input.time.size() = "<<input.time.size()<<endl;
    }
  }

  return 0;
}
